{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# IAF neurons singularity\n",
    "\n",
    "This notebook describes how NEST handles the singularities appearing in the ODE's of integrate-and-fire model neurons with alpha- or exponentially-shaped current, when the membrane and the synaptic time-constants are identical.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import sympy as sp\n",
    "sp.init_printing(use_latex=True)\n",
    "from sympy.matrices import zeros\n",
    "tau_m, tau_s, C, h = sp.symbols('tau_m, tau_s, C, h')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For alpha-shaped currents we have:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "A = sp.Matrix([[-1/tau_s,0,0],[1,-1/tau_s,0],[0,1/C,-1/tau_m]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Non-singular case ($\\tau_m\\neq \\tau_s$) \n",
    "The propagator is: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA5UAAABnCAMAAACnxJZuAAAANlBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAABHL6OuAAAAEXRSTlMAMquZdlQQ\nQN0iRInvZs27fFd4iKMAABDJSURBVHgB7V3ZoqsqDNU6nFpre/v/P3sZVKZAQAVrd3xoBUMSFoky\naaqKDhyB2wunIQpC4DAE6o84bocx/EVGXfOLtaI6fRsCd+mMVVV/mpYd3bcp+FX6NATPV7XHryoz\ncVe8fbhX1r9ax+PqdW+nZjiOHXEiBPwITOSVfnD0K+zO1U56Bp0TArkQIK+MQ/bxrqo3PSvjwCKq\nnQiQV8YByJ+T94rGlnFoEdU+BMgr4/DjY8pb28cRExUhsAsB8spd8FFhQiADAuSVEKjdezn401Gu\nIKlfqATlnYpA34xtc5mVBFRb8spTrYmEH4PAi434+/tVZuNQbckrjzEL4nImAtOTSx8vsikS15a8\n8kxrItnHIPAWu0W7zzVm43BtySsj7KJ+PZ/3+zWaPKI6v0fyEV75+Fxj5QrXlrxSt9F+bFtgX11X\nPR46GZ1/FwL9R7w68Pi036UXrE2EtuSVGnQPPl8wjFrOcjqKB+Vjmi7R7ovWf+V/+Ig2q+Xft9c6\nQlvyStWIw5PP5IET7HIe4Tb0tBVW4fU1ZxF2/jW6sts+fg8hr1TtdXuy12gmaPjYs12w7Kib12XW\nxFS1fv+sl3Z+mR6seLKHtCWvVEYrR+EqbZ+NQ9VD3VubjtKlEZDjyu4qsz1iFBzSlrxSWRA2LJnY\nsBJ6kioOdHYOAnfRl5kusjKCa0teqezoKR+ENKGjILnIWSt2ETQX2UWAa0teqQyv5a3aAysjioTO\nvhOBO5+ne15l0I9qS16pmVnbtBNNsmqAXOaULzTfruKUbHYC0Za88jKWR4r+GQTIK/9MU1NFL4MA\neeVlmooU/TMIkFf+maamil4GAfLKyzQVKfpnECCv/DNNTRW9DALklUlNRTGAkuA6hfgH2oi8Msly\nKAZQElynEP9AG5FXJlkOxQBKgusU4h9oI/LKJMuhGEBJcJ1C/ANtRF6ZZDkUAygJrlOIf6CNyCtT\nLIdiAKWgdQ7tL7QReWWK7VAMoBS0zqH9hTYir0yxHYoBlILWObS/0EbklefYDkklBPwIcK/suz//\n3Yv+/bz7QaIrhMCRCKDRfwb2EnfNZq3+/PESnzgyYKDIXAYcX5m4ZBtR9J9YW7rI99Fiq0N034sA\nRf+JbZvHRb6PFlsfovteBCj6T2zbjDSsjIWK6HYiQNF/YgF8Nd3UvmVYUk8MIIrMFQtmAboLtxFF\n/4m1j/7zYh55kwHXPDGAKDJXLJr56Xxxmq7QRhRnJNY+OjEN/eaf4A7EAKLIXLF4Zqa7dBuRV8Za\nh/zs9pN/Nt0fA6iiyFyxeGamu3QbRcQqor09woDuPJjBIB6Y/hhAFJkrs7NFs792G+GxisgruSn0\nHx7LeWTRKvoKiwFEkbminScb4bXbiKL/xBkGvzdVFYv+w6JXYjGAKDJXHKY5qa7dRhT9J8425AeY\n2AOTdWQpBlAcZmdSXbyNoqL//Pv8OxPiL5D9FjF/xkaEyqMYQF/QIogK124jNPrPf6zvRrvTERug\ny4RAUQRotqco3CSMEIhAgLwyAiQiIQSKIkBeWRRuEkYIRCBAXhkBEpEQAkURIK8sCjcJIwQiECCv\njACJSAiBoggU80o8UhJOURQZEnY6ArhF4BSnV2KLAsW8Eo+UhFNsqSCVuS4CuEXgFJesfTGvxCMl\n4RSXRJiU3owAbhE4xWbhZxYs5pV4pCSc4kygSHZ5BHCLwCnKa32ARO6VfVvgK814pCSc4oD6EosL\nIYBbBE5xoeoqVQf2+n3NP1ST+cAjJeEUmVUk9l+GAG4ROMWXVSlSnY59cbHE7nQ8UhJOoao0iXc8\nVJrOsiPQiy+NZRejC8AtAqXIbCi5QCk1rsQjJeEUa4uJ1+vWVNpJzrn0nLyNWpYRZEp5yK8WGXrk\nTeAWgVH4DMWsWVItzKL7QTH5zaqU8sqkmiPEj+eOcXDOufScvA1MygiypLTln5ZGpdMTXkOxapbC\n2Sq6GxSLn1Qlk1f2jTxYbKGq+lgHm18KXRbUAaD6555gRTnn0nPyNgApI8iW8so7bggaxRab8RuK\nXTMD3HDCLroXFJufkJ7JK8M123f15kbPSmCYcy49J2+jimUE2VKGi8Vi8RuKXTMD3HDCLroXFJuf\nkH49r6z32UbOufScvA1bKSPIkSK/LW9o8sWJgKE4NYuvhlN0JygOP65KEa9EInQgly3EgECTFkUo\nmXMuPSdvo05lBLlSAnZu6HdEAjEK5DLXwG8obs2iNXaL7gPF5cdVKeKVSPQH5LKJ2PDZtbiKzqWb\n0pJSOXkbipQRBEh58Y/LlzkQo0AuMx0DhgLULLZSQNFdoAD8mCrxXjncxTQM/5yxOpxM/vkuoOlk\nhA5VzjpDLuvUzb6Qdthcui4q9Twnb0OXMoIAKTIeqqFLMOGYB6d2MvPYTMBQgJoFq6FdBIqmgqJx\nqyqAH7se75VN3d96e0nCzuxffQUtrcq1rsc0AR7LtBCXJxasbmoe3RicYhWxQIyKUaIYAoP4xny8\nONs8REk7M5PNFDOUVFAi4Iv2Su6PPGSVcTiZt1vLn8n2MUfouA09cJFNegvGHV/zfXZVJ77KarOY\n0zIWCE/072fKY9NanZFJj5AlO1bEFt6LjKT/UoJsOauS80fL13T4xDEPTu5k5rEZZSiajna9RFq7\nHj61Sy/UcaB07+XgCNi8TE2ivZJxekDuYmY+zQ7uovf8Xzev0HPwVlc9c7RRLHJaRZck13c5/MP5\nhWL3fwERu3UsyQDciRJSwDSPmdLMzGMzuqGEFDzgWjIoqMwUr2zEg87qaMrMau6csnleeTeEBGNx\nc5hH8p0O9yrglvpo4ROgg+RvyDtORL9rkmqD6lmKtNpNMUrAaTajG0qUptuJkkFBRaV45V086KyO\npsys5s7pmzkKNK4UeiBxc3i5lj1rb23AgGU8I8HusW/hEoWGERwoooNH1DFafBFNlzoDfprNaIaS\nG79kUFCFkrxScjM7mvPYbu6c9mxcmdX8ZOw/ociYMqxkb8a8ns/7nffp449YERG8j/HKCEHxtQtQ\n+uTUqZ2HuY3K24xmKFo9ffXSSAKnntLJoHARHl5SeopXzvpCHU2sczoX3f/3UT7/YnO27Vs8V9t5\n8KwuLpL4rDuffWad4+oRHPQuJfT/sIgk3rZXelUGLiQJ0vXn5wA/QQLk43L6DzhdZ8t00uVtRjOU\nCq+Xo7DKUIV9JhQGRSuvmCLmmO6VYEcT6Zzq6uw61wDoPy/mbWK/0zjWw60eBrfn++AvdA/zLJVY\nFbWGxUFtwiLSeFte6VUZuJAmyKoRwE9QAPmunGW6QPHUjV3lYmflbUYzlMqtV4IRGIUrubBur/CF\nQDHL60AFzDHdK3XGxc8Hda/uxMzSmy2q1HzWx1m14boNbKGFvaAyT/yKVVFrWBysQVBEIm/TK70q\nAxcSBZk1AvgJAiAfkLNMFyieWJRlRXnumWYoQL3ijcAsLBfWHVQCoFjldVAC5mh4Zf+6r4dYwwBX\nVU7IXCujdeEbUatlrbgGX/67PSe2fjqPJZdFU2j9xam5kBgUEc1bKm96Jc+DVXYvRAuCK+EXZCkA\nyXHWsuDhmqiiK/5ES9EMBaqXOc4V6rPbt23+LN8oPC+sVzYqAVDM8rMgCZd8joCaGF6pF4LPjU5i\nPY5yWAfT7spdeFuLH7V6Vt55v3RdKx7B8c4H8FVoiONRNSgimvf44sf9Kf60t9BglZkq1oVoQZ5a\n2PxWMlyOO10QMMCVb9KJYVHRJS2zcMtphrILP6iwg0oAFKi8oS1ojhFeOb1vzdiL5Yre2EfAnP3B\nmALH2LybtNlOm8nCuza9TXVMerH5a3zOy6NyBn5qpprN7lRt+xa9VqBrAQ5xbOkyDYpgTtO23CgS\nebvPSlhlxnlebNooyK2LTxAux50uAGrtCsRyvBaFFVyvW2ax5q8nylBSG2plIU6g+jqoQEQzm8Al\nQQGbI+qVg3jZmu0940zehquxLmQPrhmym0mVtuHBmR1deZu7Yvt1YM0VZxv0xkp0UAe5uj09+LP0\neesr+eWGeS+UOzXLy6IHKKJqHlXNOx+JvB2v9KjMn/9Ss42CnGr5BG2So/oqjpzYjIBFxbJgt8Z5\nrsBXQhlKakOZHKNaOQBKVHlTJEthXlnP+5EnbogPx7z5J/Lcg0981knrzZy7fUje82hwubjee6TX\nM+3E47sVegy9+OcqNzKbjz17uTKycIj/h0SwO9O7E5v0xZea4nk7XulRmS1lzJhuFOTUzydoi5wh\ndb3S0YYZhrwFeyzKLQDmWGbh0qyGUiU2lMkqpnAIlJjypkSewrxy9vWq5g7pjCJ7vjjhHnxnY9pO\nesArF96yL7oIeS4jxbfo2o7S+9h+IEnA/2v+WL/LkUfbtNs/PgiLqG8fMftbpfF2vNKn8lqXjYIW\npNZ/n6AlP0XO4lIr8/STsEVF8zPNwi22Ggq7ze0xgojCQVAiyrvKI17ZLj1U8XiQN3H1fhx7tcvl\nOOd0rKSi1Mmg3MUr1bWVt7lX6C1mXnV2xjn/0FbLmLHeGdLFMYolJDreO3f6DDgDxyuXIj6Vtwpa\n+Dr/HkFJcnir7jsAi1KN7mXtkphm4RZEDMUtsD1nPyi2bMQrjd2EtfSI9f041n0zp3905ndmuSul\nfgHMXbxyLaF4d8vTUTDh0zv+QwyVuM+w7pExMeUvknqFT5EsC6ApZR/mtNVa1KvyVkErZ+vEJyhJ\nThh+SyKYBCxqbXSwgMh0SUyzcEvu19Tl6ck5XhTilcacL//GiJzylE7EP04AjisZ1ciWAPhtdXE3\nXlIeYO5Mpq4p3lLqUjp8X3pwD+a96uGW8HmDhXfU/9ROfAv9YYdX5VKCkuS8jVvkFhBci1KN7uUH\nkJhm4ZYMG4pLvyNnPyi2cMQrxYodKzPwO70wIM7AfD+O51TmO53VNK/LAZR2+RtfxxOLeeJp4pSY\nn9BCCp/09Txz5uv0lxWB/eiDFuU0OquEZVE2iWkWbqULGsp+UGz1Ea9s5oeYWH0UCwKcwfxOpc1L\npTmEDz6wA1+vW8pr2wnVI9Xhvd4LJPvjb0xKbTpDEEibWAeZgRYFm4lRfraL1WYsszBoRaKYoRwA\niq094pXVnT+aerk6xD8UII5l+XlOOn/1q+s6sbYpKe1th3N57YMhyisd3pM5QNz17SJHU8pIQmAM\nz7VF8YIsCjYTg51tM5ZZGLQiUcxQjgDFUh/zympkG3uW+JbL24nSO/0b759iAySXNPuxtdlvztW2\nE2peKRVc74nVzRzF9ftXzCwIKBmNwPOI4QNgUbCZGGrZNmOZhUErEsUM5RBQTP1Rr9TJzQ9l2U9A\nndI6Z5C63/7QtxMqr5xLqueofXsW67IWf0oWQaALToBvUcG0KNBMDLbKZmyzMMhEopChHA8KvovA\nrKzpPNYT0CTVU57Nfuw7lHxyjR/apm2ZsT5HO+f2vCv4j+ROv5sQEJ3PTSW9hXSLgs3EKLpuQXXN\nwqATiTKGkgEUdG+PWVdz0Q2/tZmlE1LrPZG3lHVM+G3SKkHJQxDIAbxpUdFqAmbhls2hryMli5Ck\nHizrh2o7ZiJubU4dYjPWe+L6NNVKvjdsrdGK0+k2BIYnuL1yG7O1lG5RayZ6ApmFW6iAoeQBRXql\nmJ1xn0puPb8ip3fmab9CrR9Xop+3Fl+omvkN5XhQZNwQ9tpQz75Lxw5zrvObwR/4Syl0lEVA7tUv\nK3OvtOyGcjwo7OMZ/Nhb8zPKD1pH+gz5f1Bmf52bttY6mQ0lIyj/AwbqsbYzxLseAAAAAElFTkSu\nQmCC\n",
      "text/latex": [
       "$$\\left[\\begin{matrix}e^{- \\frac{h}{\\tau_{s}}} & 0 & 0\\\\h e^{- \\frac{h}{\\tau_{s}}} & e^{- \\frac{h}{\\tau_{s}}} & 0\\\\- \\frac{\\tau_{m} \\tau_{s} e^{- \\frac{h}{\\tau_{s}} - \\frac{h}{\\tau_{m}}}}{C \\left(\\tau_{m}^{2} - 2 \\tau_{m} \\tau_{s} + \\tau_{s}^{2}\\right)} \\left(h \\tau_{m} e^{\\frac{h}{\\tau_{m}}} - h \\tau_{s} e^{\\frac{h}{\\tau_{m}}} + \\tau_{m} \\tau_{s} e^{\\frac{h}{\\tau_{m}}} - \\tau_{m} \\tau_{s} e^{\\frac{h}{\\tau_{s}}}\\right) & - \\frac{\\tau_{m} \\tau_{s} e^{- \\frac{h}{\\tau_{s}} - \\frac{h}{\\tau_{m}}}}{C \\left(\\tau_{m} - \\tau_{s}\\right)} \\left(e^{\\frac{h}{\\tau_{m}}} - e^{\\frac{h}{\\tau_{s}}}\\right) & e^{- \\frac{h}{\\tau_{m}}}\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡                                      -h                                     \n",
       "⎢                                      ───                                    \n",
       "⎢                                      τ_s                                    \n",
       "⎢                                     ℯ                                       \n",
       "⎢                                                                             \n",
       "⎢                                       -h                                    \n",
       "⎢                                       ───                                   \n",
       "⎢                                       τ_s                                   \n",
       "⎢                                    h⋅ℯ                                      \n",
       "⎢                                                                             \n",
       "⎢         ⎛        h            h              h              h ⎞     h     h \n",
       "⎢         ⎜       ───          ───            ───            ───⎟  - ─── - ───\n",
       "⎢         ⎜       τ_m          τ_m            τ_m            τ_s⎟    τ_s   τ_m\n",
       "⎢-τ_m⋅τ_s⋅⎝h⋅τ_m⋅ℯ    - h⋅τ_s⋅ℯ    + τ_m⋅τ_s⋅ℯ    - τ_m⋅τ_s⋅ℯ   ⎠⋅ℯ           \n",
       "⎢─────────────────────────────────────────────────────────────────────────────\n",
       "⎢                           ⎛   2                  2⎞                         \n",
       "⎣                         C⋅⎝τ_m  - 2⋅τ_m⋅τ_s + τ_s ⎠                         \n",
       "\n",
       "                                             ⎤\n",
       "                                             ⎥\n",
       "                                             ⎥\n",
       "                    0                     0  ⎥\n",
       "                                             ⎥\n",
       "                    -h                       ⎥\n",
       "                    ───                      ⎥\n",
       "                    τ_s                      ⎥\n",
       "                   ℯ                      0  ⎥\n",
       "                                             ⎥\n",
       "            ⎛  h      h ⎞     h     h        ⎥\n",
       "            ⎜ ───    ───⎟  - ─── - ───    -h ⎥\n",
       "            ⎜ τ_m    τ_s⎟    τ_s   τ_m    ───⎥\n",
       "   -τ_m⋅τ_s⋅⎝ℯ    - ℯ   ⎠⋅ℯ               τ_m⎥\n",
       "─  ────────────────────────────────────  ℯ   ⎥\n",
       "              C⋅(τ_m - τ_s)                  ⎥\n",
       "                                             ⎦"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "PA = sp.simplify(sp.exp(A*h))\n",
    "PA"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the entry in the third line and the second column $A_{32}$ would also appear in the propagator matrix in case of an exponentially shaped current"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Singular case ($\\tau_m = \\tau_s$) \n",
    "We have"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAALkAAABSBAMAAAD6AGguAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAMquZdlQQ3SJEzbvv\niWYEN0CTAAADrUlEQVRYCd2ZT0gUURzHv7ruzKq5LhR16LBKUkf/UESBITVBhtBCQZfMwSwhCrcg\njIJciOiYl46RHaJbiF2D7NAx2sAunQSPRYr9gQ4xvUVn3r95781rBgPnsr/3+73fZ96+mVm++x1A\nOvZKmQwTOw9kCJNR5+VUhpks6eUgWOeXlg3dDYIAKB/0TiamX1rq5+fyI2fwRCnMFD2vl9CjcZjX\nrP05JmvhtJjPtopzm0mv2tGLi+iYY9rF8CvwmclZ0tu7kP/FtIvhAjDm02Qs/QKtC1FnF1p/Cjl2\nSG6PpxWaiKPfn+mhE/hoqg+tP/gUO3K+E/oQzcTRaVWKHtXR9EfKRokCqZXr0RBbQndeTJPjjnRv\n0mVsRFN1/c6Qtet2ppM8YJvHmogm484+5HVXlez7mOGqxlDDlDuPou6OvAus+uFkEsc9TbQsRi2L\nyM2JSWY8AhxlhpZ03MI13cVxK87bFPTrXjfTLYXO0nH25LZrl3jahER/GAS/tR02RYneP37Fpl8/\nV6Q7aNc3WFVFOvDEql8/WaYfA1Yurjz+qO9LVpXpN4Bzczsqo8n69bNkehVoLk3ik74vWVWmkz4X\nZ1FlfoxEVGIlGEsfxzB6aiIzGidXgrH0iKMINIqE79hW9Daf/3JxI/XOsEqPdIo7U3iTiq5Xes6r\n3lR0k9I7k4puUnpJ6GolaFJ6CehqJWhUegno5GZQHEal9z/oVxtCcPoZEAo18hmzfkdWepd9dl6q\ntUNQeg+APbXs6ILSIw81/68s3dpHDEovHd2g9AY+3FxmN8oyTqn0CjNBMJv4lOJvpKkxt3u/b5pD\n67Z0H/O02RjZ0tFaNzLpBGt6bhnFQ9+OjJYoQx1Z090S8qfxLp/oK9jTyUq7nWpHRb1iWrGmN1qr\n+T43TxnqyJpOVJ7T1TLbfErNbJoPa7b0BCqvaaD6r3SotUyIRH5b0g2eHdLtjMGzS0c3eXbp6CbP\nLh3d5NkBapVH78jFMBSeJoNnB6hVXkhE4fD60OZAoBs8u4iQMNhSusGz29B5a/LKQyUomHbC2k2e\nnczVZgS6ybNDcWFwoqYlskWBbvLsCi/hvGb79bFAN3l25Vngi57IVkW6wbN7T3p3sf36WKTrZ7c0\nfFStb8abdnZ0ctENB2/a2dHJDas/BNPOju7WgcKy9gScaWdHz5GdyflaOmfaNejSuzJ1+z04BrOP\nmnYb78omPG9YzeMrHfuWfD4jjqqRaUfelXliNYOxybRLdwrWtPsLgF4+IP8NYNIAAAAASUVORK5C\nYII=\n",
      "text/latex": [
       "$$\\left[\\begin{matrix}- \\frac{1}{\\tau_{m}} & 0 & 0\\\\1 & - \\frac{1}{\\tau_{m}} & 0\\\\0 & \\frac{1}{C} & - \\frac{1}{\\tau_{m}}\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡-1           ⎤\n",
       "⎢───   0    0 ⎥\n",
       "⎢τ_m          ⎥\n",
       "⎢             ⎥\n",
       "⎢     -1      ⎥\n",
       "⎢ 1   ───   0 ⎥\n",
       "⎢     τ_m     ⎥\n",
       "⎢             ⎥\n",
       "⎢      1   -1 ⎥\n",
       "⎢ 0    ─   ───⎥\n",
       "⎣      C   τ_m⎦"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "As = sp.Matrix([[-1/tau_m,0,0],[1,-1/tau_m,0],[0,1/C,-1/tau_m]])\n",
    "As"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The propagator is"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPMAAABgCAMAAAD/2rTgAAAAP1BMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADFBd4eAAAAFHRS\nTlMAMquZdlQQQO0wRIki3e9mzbt8bDjT79EAAAX9SURBVHgB7VzZgqQqDEVFesbdGf//W28CtbKI\nQGJbc5uHri41hxyRhOKAQqSXekm3+XSLtvl0Bun+N226zadbjHJuhk8nkej/Vgk5J9pc6vJq06U+\n7lS3CrF+YjuPhqoQ1dZIKAldFNt4FAkGx28n75UzEq035Fwl1oR9uZYq0eoql89ZnK/ifZ4fiZxN\nd7D+5tX8tFJNL5vUR+1pbv8XhYtybtd7YXuUFwgMaiQLiVG4KGf7LtJ/nyfE7KlGs3G4C3BedZJs\nN6LnKA53Ac6b5txtRKkvDhfkrHopPcPLapmmcSRqEt1P1KZ/r3SbJOk2B+BCnDuMKUPvuNGKrnMO\nFh0YNl1LZT6KoND4AFyA8zBhLPUmkB4beW7aeW66ti9PMQecTLkRB+ACnOsJRmmz9xHW8bWV8AH3\npXWfhBQH8VplGpju2dYu7cEFOJtA4HNfwc8LKHUFKRUyDEHgMf25JYthOjzswYU4x9oP+OJcCcUP\njVHfxZkqV8XhApwnwzkYShXkFwnBrJblwyepxyQN1ZgkDhfgjN0VYlg5Id0RIn9GjJdTeTS8VROF\nC3AWspHzWRMhOBSA+EBVonAhzlQOXBHnh/MVW4Xep592pr+nV0T8aecrtgq9T0XtzKBOMkA6N62I\nM4M6yQBJy5lBnWSApOXMoE4yQNJyZlAnGSBJOTOokwyQDmVREsMY1EkGSFrODOokA+QRzmqdYOLn\ng0uWRre4i4LelTpLl7xr92n3iQHSOJCl0VHNQKbdA6qrszS6jmoGkopFGs4xje5r+3qF7T+7O8c1\nuj/u2ooFdBl5WxXEoNQxQL62WJZGp7YFpnhrI5DSK3Xc4l+WXtXqVUQrigsMSh0D5GsrZ+qSRlCY\nUMJgUOoYIN85H5D83LHniLLNoBubQaljgHznLOKSn8NZbaip96AhKRGVwdOVOgZIi3OGRoc3AaTl\nXoD6zKDUMUBanDM0OjMhBY0NjziDUscAaXEW6RrdqpW5vtFiLINSxwBpkf7R6Kwbgl+dGOa55l87\nhJy/fv3+12jt8vn7K2f99i7k5U/+X5/t9HX6l2/KXQd/2nn39nhOnqGoeaotPFTWzmcoaoUEPeZl\nnM9Q1DxOFx4q43yGolZI0GNexvkMRc3jdOGhIs6nKGqFBD3mRZxPUdQ8ThceCnOu+j62DZRTUaNJ\ng16UMGeY+Oy+cXafJg16UcKcYTWzelFxVGOKXphPItLtP6E0adCL4nCWtxkr9KilbedX6H3CcJYm\nDXpRHM5ifOwlUihoUJYndByVJg16URzOz+dZ1d59N3F3Q1c8oUNXPI/TpEE/isP5EbdgZ4IyuzKe\nnjz+y9pO94B+wIT/oUmDfhSHc1/Ps94ogFv/sT/jNKKzFyVrO90DOkz1cYYmDfpRHM4LbIaUWpQ0\n9atFCdxlY5eM7XQ2tA151neHM0o3ZnLbuFDXsKHO443eGpS2nc6GvqHSLCxJQLE5V7jZCZWAe5ke\nYfx+RH9mbKdzoN8AT/xic5Yw/Bo28dz5owXKYABPEekc6BNpvlVlc14hXM0rilW3suKGL09/1qeT\nttM50Pcqzv60OeMwZKh1gDKuKOjPTtjO8tKB9qNkpUEHag/F5uwY8x0ILKZJTIMZKN/HObSYBgT/\nhF3lOSjfxnlnMU1CGsxCsTh7fyQmHAz2hDcMvCq8mCYlDWahWJyDPpOfCC+muVV1KA1moRzj3Ddr\nc8vR81o3vcrb6v1qG1tMcywNZqEc4tzjwkDdy4YFx6GwxDuj3d9tY4tpjlWQhXKIM66ArDb9x4xE\nYdSSXCq9CEuPeNCWZjFNFsohzjjmHtDl+7xRlTFKsW1pFtPkoBzijO2C70CS9zlBFRyA46XeUmLr\nBcw/aDjrTBIaVBvwEZrWLLHLq6zENq9Gn9XjHYgKXwwo8dUb4dLjHg2zlDJ80d6ZEts93LRz+h2I\nB38/zHpbil7+CpUMvkmESOUlthFontMd/LbsKtHcwvU9V6dUVmKbUg/VtdXStu0KcWvU2TnvZVIl\ntlREEnAmHePQoIdBWO57PUtsE5w9dOl/jfFV4X8oDk8AAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$\\left[\\begin{matrix}e^{- \\frac{h}{\\tau_{m}}} & 0 & 0\\\\h e^{- \\frac{h}{\\tau_{m}}} & e^{- \\frac{h}{\\tau_{m}}} & 0\\\\\\frac{h^{2} e^{- \\frac{h}{\\tau_{m}}}}{2 C} & \\frac{h}{C} e^{- \\frac{h}{\\tau_{m}}} & e^{- \\frac{h}{\\tau_{m}}}\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡  -h                 ⎤\n",
       "⎢  ───                ⎥\n",
       "⎢  τ_m                ⎥\n",
       "⎢ ℯ         0      0  ⎥\n",
       "⎢                     ⎥\n",
       "⎢   -h      -h        ⎥\n",
       "⎢   ───     ───       ⎥\n",
       "⎢   τ_m     τ_m       ⎥\n",
       "⎢h⋅ℯ       ℯ       0  ⎥\n",
       "⎢                     ⎥\n",
       "⎢    -h      -h       ⎥\n",
       "⎢    ───     ───   -h ⎥\n",
       "⎢ 2  τ_m     τ_m   ───⎥\n",
       "⎢h ⋅ℯ     h⋅ℯ      τ_m⎥\n",
       "⎢───────  ──────  ℯ   ⎥\n",
       "⎣  2⋅C      C         ⎦"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "PAs = sp.simplify(sp.exp(As*h))\n",
    "PAs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numeric stability of propagator elements\n",
    "For the lines $\\tau_s\\rightarrow\\tau_m$ the entry $PA_{32}$ becomes numerically unstable, since denominator and enumerator go to zero.\n",
    "\n",
    "$1.$ We show that $PAs_{32}$ is the limit of $PA_{32}(\\tau_s)$ for $\\tau_s\\rightarrow\\tau_m$.:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAEEAAAArBAMAAADVtvhKAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAELvv3c2ZVESJZnYi\nqzKqLvLVAAABg0lEQVQ4Ea2TvUvDQBjGn7Rpm5i2huIiIkpxVOjQRV3qKDoUBCeHKkgVRYtDQQXb\nzbVbN4lLBjs5CKJLwMlZnKSQPyGCoKJWL8a8PWk+KPSFkN/zcSE5cgAyeYTNalgBr2EN6TmsIRfC\nGrFya9qr80jmxBraJDioEtdMZElwUCbeBOZJcLCkG39qEckOFxB2FM1h4Q1yTiCfQJqMGI4YshBt\nnFNAoGhiwhGpIobrDQoI4kbq0hERDfI2+YOAgywbFd801iCe2t8zRm933QXplfyCy927tIx7VyWk\nuovc/a4i7JGUVehH+ukZGTY8tG5UMmLAcTFaOWSG++EWvihmkAFS5giavIf3fwoQUUWhwpuzgMLr\ncVzjyeCdKwj0tbzfZWXnoiv6pbGZtkZ75bVYmQOq7PIdYYrtU/zTNwdi9sGWcgGNfXtT0iX/RvLF\nzgTTvyF/+GdOInqeSH5V1LKVx6mkkvjbOCHdC85OlHoDcpL2b8H+uYCpaUisB+TsJTfaW0GFH7To\nVTvpgivnAAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$$\\frac{h}{C} e^{- \\frac{h}{\\tau_{m}}}$$"
      ],
      "text/plain": [
       "   -h \n",
       "   ───\n",
       "   τ_m\n",
       "h⋅ℯ   \n",
       "──────\n",
       "  C   "
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "PA_32 = PA.row(2).col(1)[0]\n",
    "sp.limit(PA_32, tau_s, tau_m)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$2.$ The Taylor-series up to the second order of the function $PA_{32}(\\tau_s)$ is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhoAAAA4BAMAAABJbEwBAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAELvv3c2ZVESJZnYi\nqzKqLvLVAAAIuUlEQVRoBeVZfYwkRRV/PV89PdPTM3JAQtDseDFKLnysEYMiklFRQyQwRhaIEhgI\neKfewWjMGdTIBMGvRG/5gzNIgFbixnCGrBfuxN0jNB8xQo7sEILGXO7oQMREIbuXGF3u8NZXVV1f\n3dUz23OzxGVfslWvfu/3ftVV013dVQswPqu1x6e1/pXszvofw/hGUOqOT2v9KxXmtq//QYw0gt3v\nbyTy6uFcAtsQgBMYlsw74OwNMfjEIAuB91YCnIPWxlw58j13OTEb2+FQ8vFJsN6RgH30HTksOOX8\nTOOK6IVZLau+IkzD11/jqmyXzOi7oyR7M9p7syn8X7P/k+3qKL3SyZa0btjlbCsAo/8YfrZuBpjp\nQp1Wdrp39bfexqej4me6wpHIVp+lFbfuMX41vWRWpfT8ysoSgHXf/oCT/vahSxrcH3e9b9yCJr27\nGDhxLRw2hXeaQACFfmfgTnNS90fcG3td/uDYJQ2C1T4Fbw9hsyEKW00ggKS7F8DpoSDdDZWtD2+7\nTQIicrLOppE0TTupgVfyRRq9GeBjJtrnZgITLOm1C+d/KRltsPfBBba4WWTkZL0/jyJg3EkNFFqg\nT/ql4J4w0U5UfRMs6XV2bzGS1QQ4YrWcrinnpDAv2zsv6su4kxp4HTUyHGsZnEkrySs3c0ESVen1\n6Vi8ZXcKdgw7+WaxM4qGcSc1UIjOemUJ8tPKHc8zqr5xZAo9j7Mxy+lYW/1qr3ZAAcbjLgSj6WTe\nSX0J+8GDinov/isjXgpqvzNchkKvtuH0wEAZM3TTiHqxndRwlYeQkvPB+cpwKmeo9P1zPoeT9T+S\nUCrySDJSvuGMK66hMFviPT/JSUMqIYnwnVQaK4Ev9BIQA8SGlHxkjWRuM0NaLUyQHy+f9dpKD+Fo\niX8qwRgAHMFY9p0UefJ1+wbZkjZ0bKRWMZNIO95Hbfb7DXiCwLhQEfs8LVdZFJGXfSdlONxM6+/d\naYEU/A8pOIUTYpxtBVHaDusjAHUyEaUWgarTpEyxhByyR9hJVdm8p3Siwd/RWsMbFw2iJMQ2RWx7\nljnlo6UmQJF8aRTQAcgHpEyxhJw1CWwnlZIAB0wB+y0TasQSPRpZAiy3hGtwEmLFLmPx2Sh0bARy\n5Meqd0joTlKkWUIOrk6jRrhxD+ZpBz1nzIv/DnlfOP8TumCyRz0eazkdCVRxTT4um+hFYjLA6Xw2\ndoWEn5vEok5vF/UF8QTqUZBwiPFrkwFBL/QxLnFKx6LAHD1gHeNhrMuXwbO8aZd73I1q3mMMTmvi\ni1jYva8fORiIFnEiMRmotFmcz8atDdLOd7CY8LGALaRg5jzw/OsHeYPWkZwS+DWPb0JhBeew3SRe\nPHCch7F+umvdIppOA2a+PfPDBwSQcTYKXZHpdaEpGsxhYkrAbrEAn43zaHOxh9WCT/zPkoLZ3yHH\n3aiOrk0JfFNlKLiA95LpjgfeBLzrqC3Bi3vm6C9CM4r4A7bz3dtEOp8Nzh9QL2FSvgdw2iXErgPw\n+kKGOVxMBDySg8Zn43p6JTcTaLFLSnzDWFcSuU+HAG8QRDEuJwOLShRdmVAiGmhXLFOGDJDmvynG\niv8qPpyCn+zhqfBzht0zNfWpqamrVMIQn8yGMMcXLjqKmAjQ2XCmpq788NRUBzkvh1i4dI1f9NEl\nsyHtF9LV5EAEYrMhcJF4V4+6ekCdjTcFlTkF2Kn8K1HMf4yV0lSeFFyzQqh89Z/7q3vDiM3FMADO\nPP745RaL8HvjwT62T+sRMPGkAFwOMPPIfGf+XBJH43IkAPf9HovYK0jghE3Miz7mdCV1Bj6K3ziM\ny8r3wGNwKOCI6JED8Vo/a1JXUVjAL8pdfv6FKhkjMS6GAdjmNTHcJrB8UiaO4a1xFoUmZkm1hfqs\ncJcAfnMMPt54PAK5HAmUeq8gKlZRypB4lADFQAQUJexU2D6wtotGwuE9JgIREDtrcvoK8VGAJx/F\nnXJpOgK5GAZg89cbOIYOi/B7o7Lylz/tCCjGTlHuZ3FaVlvgntmCZeCbTS6HAch94F1I+jKnb2qi\nJ3EO/4A5upKrfm9Uv/YwJxtq3qMhRKHYWRO/9WnsEJaH4FUQjw8XI4E3VkL86uyhh8ZnA+5Y+UxI\nEcjTeVVv/dosLmqz3qR40XA5EvBeJgd517BcfEgn0ZM4h29hjq7kLRP0zHMO+89yXlrNe2TxP954\nLnGUzPhZ02UxoYthLyx4EaiI2QG5pKejgJgNmZ1rE78QSIR4+bDadJdCBipy8F2cd3DJHDBDcWoU\nj3ytkkr0eKh6IcBO/BtsP1HDVh8e7OEdqGXqZ00zKh+vbwmv8ldOBCpihbAU4tIXBcp+5MjKoSOr\ndCRCvAX8wSvNaQYqcriQ4u/u8DmQh3MU1zVYSyqRnqz3NfC5Xf2GhWjUGpCfjGfqZ025gHUWleU+\nXAvP361htFGZewyvoZ0McCSa5LN5m9U/xYXQ/Steedz27Onh20iglh+5FBew4kilXD/aKJYnlfhw\nN9+C2gm2xZSZu7U8t6k1BzdiU6eRo93DUxo4pCGnrjKEqYbrPsCtXUS8jgoP9XGrjbOhZ8bv5eeG\nqkjCb6Wb9M6hkOsnI2lIJUiLDMQXA3D/RRhWOJBnCBaOiky2Fc1+1mRQNUEPmcC1wO7BBef4aMK7\npkUm3YqOcNa0yo4nyM37dhi+Ewrk9TyCfVJksq3o4LOmEToQKdW2cNfUsXHtzC+RLqys/ZSaSqbX\nz5qeiX9xJvbI5IKP9wadje9l1cAtgsx0UGcNTd90rFlHpBv2pdHJ2IfdgYMyk2xF19Bys2soLqTd\ni9B1yTbWaQhwdc5rAAdkJtmKrqFZl6+huJA+tUvc232wbxDY6hz30nt3tGQmbkX184vVqaya9Uyw\nauroxC001brx8LasGkU8B2yByMR9kn5+kVVvCN8NhxDGELZ6YxDhEvr5BUc3bK2dX2zYWeAD184v\nOLhRa/38YqPOAh936vkFJ6yL+n+kem1jDpDuygAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$$\\frac{h}{C} e^{- \\frac{h}{\\tau_{m}}} + \\frac{h^{2} e^{- \\frac{h}{\\tau_{m}}}}{2 C \\tau_{m}^{2}} \\left(- \\tau_{m} + \\tau_{s}\\right) + \\mathcal{O}\\left(\\left(- \\tau_{m} + \\tau_{s}\\right)^{2}; \\tau_{s}\\rightarrow\\tau_{m}\\right)$$"
      ],
      "text/plain": [
       "   -h                     -h                               \n",
       "   ───                    ───                              \n",
       "   τ_m    2               τ_m                              \n",
       "h⋅ℯ      h ⋅(-τ_m + τ_s)⋅ℯ       ⎛            2           ⎞\n",
       "────── + ──────────────────── + O⎝(-τ_m + τ_s) ; τ_s → τ_m⎠\n",
       "  C                   2                                    \n",
       "               2⋅C⋅τ_m                                     "
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "PA_32_series = PA_32.series(x=tau_s,x0=tau_m,n=2)\n",
    "PA_32_series "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Therefore we have \n",
    " \n",
    "$T(PA_{32}(\\tau_s,\\tau_m))=PAs_{32}+PA_{32}^{lin}+O(2)$ where $PA_{32}^{lin}=h^2(-\\tau_m + \\tau_s)*exp(-h/\\tau_m)/(2C\\tau_m^2)$\n",
    " \n",
    "$3.$ We define\n",
    "\n",
    "$dev:=|PA_{32}-PAs_{32}|$\n",
    " \n",
    "We also define $PA_{32}^{real}$ which is the correct value of P32 without misscalculation (instability).\n",
    " \n",
    "In the following we assume $0<|\\tau_s-\\tau_m|<0.1$. We consider two different cases\n",
    " \n",
    "a) When $dev \\geq 2|PA_{32}^{lin}|$ we do not trust the numeric evaluation of $PA_{32}$, since it strongly deviates from the first order correction. In this case the error we make is\n",
    " \n",
    "$|PAs_{32}-PA_{32}^{real}|\\approx |P_{32}^{lin}|$\n",
    "        \n",
    "b) When $dev \\le |2PA_{32}^{lin}|$ we trust the numeric evaluation of $PA_{32}$. In this case the maximal error occurs when $dev\\approx 2 PA_{32}^{lin}$ due to numeric instabilities. The order of the error is again\n",
    "\n",
    "$|PAs_{32}-PA_{32}^{real}|\\approx |P_{32}^{lin}|$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The entry $A_{31}$ is numerically unstable, too and we treat it analogously."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Tests and examples\n",
    "We will now show that the stability criterion explained above leads to a reasonable behavior for $\\tau_s\\rightarrow\\tau_m$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import nest\n",
    "import numpy as np\n",
    "import pylab as pl"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Neuron, simulation and plotting parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "taum = 10.\n",
    "C_m = 250.\n",
    "# array of distances between tau_m and tau_ex\n",
    "epsilon_array = np.hstack(([0.],10.**(np.arange(-6.,1.,1.))))[::-1]\n",
    "dt = 0.1\n",
    "fig = pl.figure(1)\n",
    "NUM_COLORS = len(epsilon_array)\n",
    "cmap = pl.get_cmap('gist_ncar')\n",
    "maxVs = []"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Loop through epsilon array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x598b510>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "for i,epsilon in enumerate(epsilon_array):\n",
    "    nest.ResetKernel() # reset simulation kernel \n",
    "    nest.SetKernelStatus({'resolution':dt})\n",
    "\n",
    "    # Current based alpha neuron \n",
    "    neuron = nest.Create('iaf_psc_alpha') \n",
    "    nest.SetStatus(neuron,{'C_m':C_m,'tau_m':taum,'t_ref':0.,'V_reset':-70.,'V_th':1e32,\n",
    "                           'tau_syn_ex':taum+epsilon,'tau_syn_in':taum+epsilon,'I_e':0.})\n",
    "   \n",
    "    # create a spike generator\n",
    "    spikegenerator_ex=nest.Create('spike_generator')\n",
    "    nest.SetStatus(spikegenerator_ex,{'spike_times': [50.]})\n",
    "    \n",
    "    # create a voltmeter\n",
    "    vm = nest.Create('voltmeter',params={'interval':dt})\n",
    "\n",
    "    ## connect spike generator and voltmeter to the neuron\n",
    "\n",
    "    nest.Connect(spikegenerator_ex, neuron,'all_to_all',{'weight':100.})\n",
    "    nest.Connect(vm, neuron)\n",
    "\n",
    "    # run simulation for 200ms\n",
    "    nest.Simulate(200.) \n",
    "\n",
    "    # read out recording time and voltage from voltmeter\n",
    "    times=nest.GetStatus(vm)[0]['events']['times']\n",
    "    voltage=nest.GetStatus(vm)[0]['events']['V_m']\n",
    "    \n",
    "    # store maximum value of voltage trace in array\n",
    "    maxVs.append(np.max(voltage))\n",
    "\n",
    "    # plot voltage trace\n",
    "    if epsilon == 0.:\n",
    "        pl.plot(times,voltage,'--',color='black',label='singular')\n",
    "    else:\n",
    "        pl.plot(times,voltage,color = cmap(1.*i/NUM_COLORS),label=str(epsilon))\n",
    "\n",
    "pl.legend()\n",
    "pl.xlabel('time t (ms)')\n",
    "pl.ylabel('voltage V (mV)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Show maximum values of voltage traces"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x6799350>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fig = pl.figure(2)\n",
    "pl.semilogx(epsilon_array,maxVs,color='red',label='maxV')\n",
    "#show singular solution as horizontal line\n",
    "pl.semilogx(epsilon_array,np.ones(len(epsilon_array))*maxVs[-1],color='black',label='singular')\n",
    "pl.xlabel('epsilon')\n",
    "pl.ylabel('max(voltage V) (mV)')\n",
    "pl.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "pl.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The maximum of the voltage traces show that the non-singular case nicely converges to the singular one and no numeric instabilities occur. "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
